Dielectric breakdown in spin polarized Mott insulator 
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Nonlinear response of a Mott insulator to external electric field, corresponding to dielectric breakdown phe- 
nomenon, is studied within of a one-dimensional half-filled Hubbard model. It is shown that in the limit of nearly 
spin polarized insulator the decay rate of the ground state into excited holon-doublon pairs can be evaluated nu- 
merically as well to high accuracy analytically. Results show that the threshold field depends on the charge gap 
as F t h oc A 3 / 2 . Numerical results on small systems indicate on the persistence of a similar mechanism for the 
breakdown for decreasing magnetization down to unpolarised system. 
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The nonlinear response to external fields and more general 
nonequilibrium properties of strongly correlated electrons and 
Mott insulators in particular fl]] are getting more attention in 
recent years, also in connection with powerful novel experi- 
mental techniques, e.g. the pump-probe experiments on Mott 
insulators J2], as well as novel systems, the prominent ex- 
ample being the driven ultracold atoms within the insulating 
phase [3]. In this connection, one of the basic phenomena 
to be understood is the dielectric breakdown in Mott insu- 
lators, studied experimentally in effectively one-dimensional 
(ID) systems more than a decade ago fl. The concept of 
Landau-Zener (LZ) single-electron tunneling HE] as a stan- 
dard approach to dielectric breakdown of band insulators |7] 
is not straightforward to generalize to correlated electrons [8- 
rToh . Theoretical efforts have been so far restricted to the pro- 
totype Hubbard model at half-filling. In ID numerical ap- 
proaches have given some support to analytical approxima- 
tions for the most interesting quantity being the threshold field 
F t h and its dependence on the charge gap A [9], typically re- 
vealing a LZ type dependence F t h oc A 2 . Different depen- 
dence is found numerically within the dynamical-mean-field- 
theory approach 11 10 as relevant for high dimensions D 3> 1. 



In this Letter we approach the problem of a dielectric break- 
down from a partially spin polarized Mott insulator. We use 
the fact that the ground state (g.s.) of the ID Hubbard model 
is insulating at any spin polarization with the charge gap mod- 
estly dependent on the magnetization to. In particular, a sin- 
gle spin excitation in fully polarized system to ~ 1/2, i.e. 
AS* = 1 state, can be studied exactly numerically as well as 
to high accuracy analytically. The relevant mechanism for the 
decay of the g.s. under constant external field F is the cre- 
ation of holon-doublon (HD) pairs. We show that due to the 
dispersion-less g.s. the similarity to the LZ tunneling is only 
partial and leads to a different scaling F t h oc A 3 / 2 . Furtheron 
we study numerically on small systems also the model with 
AS > 1, m < 1/2 in a finite field F. Results indicate that 
the decay mechanism remains qualitatively and even quantita- 
tively similar at polarizations m < 1/2, in particular for larger 
A whereby the most interesting case is clearly the unpolarized 
m = system. 



In the following we study the prototype ID Hubbard model, 

H = -tJ2(^4+i,a^ + H.c) n n n iU (1) 

■ia i 

with periodic boundary conditions (p.b.c.) where cJ CT , Ci„ are 
creation (annihilation) operators for electrons at site i and spin 
(7 =t, 4- The action of an external electric field F is induced 
via the Peierls phase <f> (vector potential) and its time depen- 
dence, i.e. qb(r) = eoF(T)ao/h. Furtheron we use units 
h = erj = clq = 1, as well as we put t = 1 defining the 
unit of energy. In such a model we investigate finite systems 
of length L and at half-filling N u + Nd = L but in general 
at finite total spin, S z — (N u — Nd)/2 and magnetization 
to = S z /L. 

Let us first consider the problem of a single overturned spin, 
i.e. AS Z = L/2 — S z = 1. Here, basis wavefunctions \(pj m ) 
correspond to an empty site (holon) at site j and a doubly oc- 
cupied site (doublon) at site to. Taking into account the trans- 
lational symmetry of the model (Q]) with p.b.c. (even with time 
dependent 4>(t)) at given (total) momentum q = 2irm q /L 
the relevant basis is | = (1/vT) X] ■ e iq3 : \ipj t j + i),l € 
[0, L — 1]. At fixed qb adiabatic eigenfunctions can be then 
searched in the form \tp) = J^j ^jl^q) leading to the eigen- 
value equation, 

1 _ 1 ^ 1 

~ U ~ L^f E -U + 2(cos(g' -(/))+ cos(q' -<j>-q))' 

(2) 

In the limit L —> oo the g.s. energy E representing the 
holon-doublon (HD) bound state can be expressed explicitly 
as E Q = U - {U 2 + 16 cos 2 (g/2) j 1 / 2 . We note that (in 
spite of the g-dependence) g.s. states for all q are noncon- 
ducting since from Eq. (|2] ) it follows that the charge stiff- 
ness Vq oc d 2 E /d(j) 2 — ¥ for L — >• oo. On the other 
hand, excited states form a continuum with lower edge at 
E x = E/-4cos(g/2). 

Since qb(r) conserves total q we furtheron consider only 
solutions within the q = subspace representing the ab- 
solute g.s. wavefunction |0) with <f- = Ae~ K '^'e % ^ and 
A = Vtanh n. Here, the charge gap A = E\ — Eq and 
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the related g.s. localization parameter k are given by 



A = -4 + yJ\J 2 + 16 = 4(cosh k - 1). 



(3) 



When we consider the time-dependent <p(r) we have to deal 
at finite L with adiabatic states E n ((f>) as, e.g., shown in 
Fig. 1 for finite L. At finite L ^> 1/k Eq is essentially <p- 
independent but the same holds as well for lowest excited 
states E n , n > 1 which makes an usual application of two- 
level LZ approach not straightforward to apply. If we sup- 
pose that due to field F > the transition probability be- 
tween neighbouring states is high (neglecting finite size gaps 
between them) excited states are well represented by 'free' 
HD pair states with 



4 = ^= e 

3 Vl 



e k = U — 4cos(0 — k), k 



2n 



-m k . (4) 



As shown further relevant transitions due to time-dependent 
(/>(t) happen to effective states with \irik\ 3> 1 since the 
g.s. |0) is well localized. 




Figure 1. Energy levels E n (in units of t) vs. phase (f> for holon- 
doublon pair states in the system with L = 21 sites and U — 4. 
Thick line represents the g.s. and the effective HD pair state disper- 



Let us now consider the decay of the g.s. |0) after switch- 
ing constant field F(r > 0) = F,<f> = Ft. We present an 
analysis for the initial decay where most weight is still within 
the g.s., i.e. |ao(r)| ^> |a n ^n(T)|. In sucn case the excited 
state amplitude time-dependence a„(r) is given by 



a n (r) = -F [ dr'$„(r')exp(z f w n (r")rfr"), 
Jo Jo 



(5) 



where $„ = (n\d/d(j>\0) and w„(t) = E n {4>) — Eq. 

Analytically progress can be made by using effective HD 
states \k) as approximate excited states with u>k(0 — e fc — 
Eq = 4(cosh«; — cos 5), £ = Ft — k. By using the relation 

(fc|0)w fc = (k\H + U- H\0) = U(k\n oi \0) = UA/VZ, 

(6) 

where Hq denotes only kinetic term in Eq. (Q]), one can express 
in Eq. (O as 



Here, we can already realize some essential differences to the 
usual concept of of LZ tunneling, i.e., <£>£ and Eq. © do 
not favor transitions to lowest lying excited state but rather 
to k ~ k/v3, hence the reduction to a two-level problem is 
not appropriate. 

The rate of afc(r) following from Eqs. (0, (0 is not steady. 
Since we are interested in low F we average it over the Bloch 
period tq = 2ir/Fto get a = a j, (tb ) which is approximately 
the same for majority of k (fixing here k = tt), 

4 =-$ /X=ao)' exp (i£^' (r) ) <8) 

after per partes integration of Eq. ([8]) and neglecting the first 
fast oscillating part, smaller also due to an additional prefactor 
F. Final simplification for small F can be made by replacing 
coshK — cos£ ~ £ 2 /2 + A/4 and consequently extending 
integrations in Eq. (|9]l to £ = ±oo. This leads to an analytical 
expression for the decay rate r, defined by |ao| 2 ~ exp(— Tr) 
where V = L\a\ 2 /tb, 



r 



A 3 / 2 B(A) 2 / y/2A^ 
-Ki 1 



3nF 



3F 



B(A) 
V8 



■ exp 



(2A) 



3/2 



3F 
(10) 

where ^1/3(0;) is the modified Bessel function and £>(A) = 
A(A + 8) 3 / 2 /(A + 4), and the last exponential approxima- 
tion is valid for small enough T. The main conclusion of 
the analysis is that F in Eq. (TTOb depends on A 3 / 2 /F un- 
like usual LZ theory applications [8, 9] yielding A 2 /F. As 
the threshold field is usually defined with the expression T oc 
exp(-7rF tfe /F),Eq.(0 directly leads F th = (2A) 3 / 2 /(3tt). 

It is straightforward to verify the validity of approxima- 
tions for Nd = 1 via a direct numerical solution of the time- 
dependent Schrodinger equation (TDSE) with cj> = Ft within 
the full basis at q = and finite but large L > 100. Time de- 
pendence of the g.s. weight |ao(r)| 2 is presented in Fig. 2 for 
typical case {7 = 4 and different fields F = 0.2 — 0.5. Results 
for the case of an instantaneous switching F(t > 0) = F 
(shown for F = 0.5) reveal some oscillations (with the fre- 
quency proportional to the gap A) but otherwise clear expo- 
nential decay with well defined T. In order to minimize the 
fast-switching effect we use in Fig. 2 and furtheron mostly 
smooth transient 111 ill , i.e., field increases as F(r < 0) = 
Fexp(3r/r B ) to its final value F(t > 0) = F. 

In Fig. 3 we compare results for V as obtained via three dif- 
ferent methods: a) direct numerical solution of TDSE, b) ana- 
lytical approximation with an average decay rate into free HD 
states, numerically integrating Eq. (©, and c) the explicit ex- 
pression ( [Tol l where additional simplification of the parabolic 
dispersion of excited states is used. The agreement between 
different methods is satisfactory essentially within the whole 
regime of small T and deviations between analytical and nu- 
merical results become visible only for large T ~ 0.1. More- 
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Figure 2. (Color online) a) Ground state weight In ao| 2 vs. time 
t/t b for U = 4 and different fields F = 0.2 - 0.5. For F = 0.5 
the comparison of results for smoothly and instantaneously switched 
F(t) is presented while for F < 0.5 only smooth switching is used. 

over, results confirm the expected variation In T oc 1/F es- 
sentially in the whole investigated range of F. 




Figure 3. (Color online) Ground state decay rate T (log scale) vs. 
1/F for U — 4 as evaluated by direct numerical solution of TDSE 
(full line), decay into free HD states, numerically integrating Eq. ([8} 
(dotted line), and analytical expression, Eq. J 10b (dashed line). 

One can assume that a similar mechanism of the dielectric 
breakdown via the decay into free HD pairs remains valid at 
finite deviations N4 > 1 and m < 1/2. In order to test this 
scenario we perform the numerical solution of TDSE for the 
model, Eq. (Q~|i, with the finite field F(t). Calculation for all 
S z sectors covering the whole regime < m < 1 /2 are per- 
formed on finite Hubbard chains with up to L = 16 sites using 
the Lanczos procedure both for the determination of the initial 
g.s. wavefunction |0) as well as for the time integration of the 
TDSE 1 14] within the full basis for given quantum numbers 
Na, N u , q = reaching up to N st ~ 10 7 basis states. We use 
everywhere smooth transient for the field F(t). Since the de- 
cay rate of the g.s. weight |ao| 2 is expected to scale with the 
number of overturned spins Nd the relevant quantity to follow 
and compare is (l/Nj) In |ao| 2 (r). 

In Figs. 4,5 we present numerical results for time depen- 
dence of normalized g.s. weight In \ ao\ 2 /Nd as obtained via a 
direct solution of the TDSE for L = 16 with the whole range 



of magnetization 1/2 > m > (relevant 1 < Nd < L/2) for 
two cases of U = 4, 10, respectively, and the span of appro- 
priate fields F. Examples are chosen such to represent charge 
gap (for a single HD pair) being small A ~ 1.3 < W and 
large A ~ 6.5 > W, respectively, relative to the noninteract- 
ing bandwidth W = 4. 




t/t b 

Figure 4. (Color online) Normalized g.s. weight (l/Nd) In |ao| 2 vs. 
time t/tb for U = 4 and fields F = 0.3, 0.6, 0.8, for various spin 
states 1 < N d < L/2. 
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Figure 5. (Color online) The same as in Fig. 4 for U — 10 and 
F = 1.6,2.2,2.6. 

The main conclusion following from Figs. 4,5 is that the 
g.s. weight I ao 1 2 indeed decays proportional to Nd confirming 
the basic mechanism of the field-induced creation of (nearly 
independent) HD pairs. The decay rate Y defined as |a | 2 oc 
exp(— YNdT) is only moderately dependent on Nd and m. 
Results confirm that Y is essentially independent of Nd in well 
polarized systems with m > 1/4, which is compatible with 
independent decay into low concentration of HD pairs. For 
larger U — 10 in Fig. 5 the invariance of Y extends even to 
unpolarized situation m = (Nd/ L =1/2) for intermediate 
fields F > 2.2. 

There are some visible deviations at m < 1/4 for weak- 
est fields both in Fig. 5 for F — 1.6 and even more for 
smaller U = 4 and F = 0.3 in Fig. 4, indicating on larger 
r and correspondingly faster decay of unpolarized g.s. with 
m = relative to nearly saturated m ~ 1/2. Part of 
this enhancement of Y can be attributed to the dependence 
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Figure 6. (Color online) Threshold field F t h vs. charge gap A for dif- 
ferent magnetizations to ~ 1/2 (given by Nd = 1) and m = 1/4, 
as obtained numerically for L = 16. Full curve (HD) represent the 
analytical approximation, Eq. d 10b . while the dashed curve is the LZ 
approach result from Ref. 1 9]. 



of the charge gap on the magnetization A (to). The ther- 
modynamic (L — > oo) value Ao = A(m = 0) is known 
via the Bethe Ansatz solution given by the equation Ao = 
(16/t/) f^dxVx* - l/smh(2-!tx/U) HEl]- Values for 
A(m ~ 1/2) as given by Eq. are somewhat larger than 
Ao with the relative difference becoming more pronounced 
for U < 4. Still taking into account actual A(m) some en- 
hancement seems to remain at to ~ at least for weaker 
fields F and smaller U. This could indicate that the decay 
into HD pairs are not independent processes but correlations 
due to finite concentration of Nd/L enhance decay. 

Finally let us consider the threshold field for the decay F t h 
as defined again by T oc exp(— irFth/F). We present results 
in Fig. 6 for F t h as function of the gap A. To extract F t h 
vs. A we use numerical data for T(F) obtained from numer- 
ical |ao| 2 (r) as, e.g., shown in Figs. 2,4,5. For the reference 
charge gap A(to) we use for m ~ 1/2 and m = 1/4 Eq. (0, 
while for m — we use exact Ao. Some deviation between 
to ~ 1/2 and m = 1/4 results can be still attributed to ac- 
tually slightly smaller gap for the latter magnetization. For 
comparison we plot also the analytical result emerging from 
Eq. ([Tol l. Fth oc A 3 / 2 , as well as the dependence following 
from the LZ approach H with F th = A 2 /8. From Fig. 6 
we conclude that the general trend F t h(A) is quite well rep- 
resented by the single HD pair result which deviates signif- 
icantly from the LZ dependence at least for larger A > 6. 
At the same time, we should note that our numerical results 
in the range 1 < A < 2.1 agree also well with data analyz- 
ing numerically the g.s. decay using the t-DMRG method (at 
to = 0) for the same model but bigger L ~ 50 |9[]. 

In conclusion, we have presented an analysis of the di- 
electric breakdown within the Mott-Hubbard insulator start- 
ing from a spin polarized ground state. Such an approach has 
clearly an advantage since the problem can be solved up to de- 
sired accuracy numerically but as well captured analytically. 
As such the situation can serve at least as well controlled test 
for more demanding situations of an arbitrary magnetization, 



in particular of an unpolarized g.s. |T(J 11 ]. 

The case of a nearly polarized state Nd = 1 describes the 
mechanism of the field-induced decay of the g.s. into sin- 
gle HD pair. Here one can follow differences to usual LZ- 
type approaches: a) the g.s. is localized and dispersionless 
within the insulator, b) the transition is not between two iso- 
lated levels but rather to a continuum, moreover it follows 
from Eqs. (HJ,© that matrix elements do not favor transi- 
tions to lowest excited states, c) instead of exact excites states, 
one can well use effective free HD states, d) dispersion of ef- 
fective HD states is unlike in LZ applications not hyperbolic, 
e.g., ojk oc (fc 2 + k 2 ) 1 / 2 but rather parabolic lu^ = k 2 + k 2 
which is presumably the main origin for qualitatively differ- 
ent behavior of the threshold field Fth °c A 3 / 2 which is a final 
manifestation of the distinction to usual LZ applications. On 
the other hand there are some similarities. In particular the 
analytical expression for the average transition rate, Eq. (O, 
where matrix element is integrated out, appears analogous to 
two-level problem and ready for phase-integral transformation 
into imaginary plane as used originally by Landau [5] and then 
generalized II171 ll 811 and applied as well to breakdown prob- 
lem OlOi 1 1311 . Still it is straightforward to verify that for the 
levels under consideration Uk do not satisfy criteria for its ap- 
plication, but the analogy rather emerges through the applica- 
tion of the steepest descent approximation to Eqs. ©. 

The picture of the decay of the driven Mott insulator into 
HD pairs remains attractive for magnetization approaching 
the unpolarized g.s. There seem to be two characteristic 
length scales controlling the mechanism, the HD pair localiza- 
tion length ( = 1/k and the Stark (Bloch) localization scale 
Ls = 8/F. Our results indicate that for larger A (small Q and 
well localized HD pairs the mechanism of decay into nearly 
independent HD pairs remains at least qualitatively valid. On 
the other hand, we find indications that for smaller A and 
weaker F (larger Ls), the decay is enhanced, i.e., pointing 
into the direction of more collective driven excitations favored 
also in the interpretation of experiments [4]. It should be as 
well pointed out that the phenomenon of HD pair generation 
is not particularly specific to ID systems discussed here but 
can generalised to higher dimensional Mott insulators as well. 
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